Matrix product decomposition and classical simulation of quantum dynamics 

in the presence of a symmetry 
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We propose a refined matrix product state representation for many-body quantum states that 
are invariant under SU(2) transformations, and indicate how to extend the time-evolving block 
decimation (TEBD) algorithm in order to simulate time evolution in an SU(2) invariant system. 
The resulting algorithm is tested in a critical quantum spin chain and shown to be significantly 
more efficient than the standard TEBD. 

PACS numbers: 



i 

c/2 



I 

S3 

o 
o 



> 

o 
o 



c : 
o 

o ■ 



X 



Quantum many-body systems are described by a large 
Hilbert space, one whose dimension grows exponentially 
with the system's size. This makes the numerical study 
of generic quantum many-body phenomena computa- 
tionally hard. However, quantum systems are governed 
by Hamiltonians made of local interactions, that is, by 
highly non-generic operators. As a result, physically rel- 
evant states are atypical vectors in the Hilbert space and 
may sometimes be described efficiently. Systems in one 
spatial dimension offer a prominent example. Here the 
geometry of local interactions induces an anomalously 
small amount of bipartite correlations and an efficient 
representation is often possible in terms of a trial wave 
function known as matrix product state (MPS) 0, Q. 
This, in turn, underlies the success of the density ma- 
trix renormalization group (DMRG) an algorithm to 
compute ground states, and of several recent extensions 
@j Hi @j) including the time-evolving block decimation 
(TEBD) algorithm to simulate time evolution 

Symmetries, of fundamental importance in Physics, re- 
quire a special treatment in numerical studies. Unless ex- 
plicitly preserved at the algorithmic level, they are bound 
to be destroyed by the accumulation of small errors, in 
which case significant features of the system might be 
concealed. On the other hand, when properly handled, 
the presence of a symmetry can be exploited to reduce 
simulation costs. Whereas the latter has long been re- 
alised in the context of DMRG [!, 0] , the subject remains 
mostly unexplored for the TEBD algorithm [8j]. 

In this letter we undertake the study of how to enhance 
the MPS representation and the TEBD algorithm in sys- 
tems that are invariant under the action of a Lie group 
Q . We present an explicit theoretical construction of a 
refined MPS representation with built-in symmetry, and 
put forward a significantly faster TEBD algorithm that 
both preserves and exploits the symmetry. For simplic- 
ity and concreteness, we analyse the smallest non-abelian 
case, the SU(2) group, which is extremely relevant in the 
context of isotropic quantum spin systems. The analysis 
of the SU(2) group already contains the major ingredi- 
ents of a generic group Q — in contrast with the case of 



an abelian U(l) symmetry [8(. In addition, it can be 
cast in the language of spin operators, more familiar to 
physicists than group representation theory. As a test, 
we have computed the ground state of the spin- 1/2 an- 
tiferromagnetic heisenberg chain, obtaining remarkably 
precise two-point correlators both for short and long dis- 
tances. 

In preparation to describe the SU(2) MPS, we start 
by introducing a convenient vector basis and discuss a 
bipartite decomposition of states invariant under SU(2). 

Total spin basis.— Let V be a vector space on which 
SU(2) acts unitarily by means of transformations e lv ' s , 
where matrices S x , S y and S z close the Lie algebra su(2), 



namely [S a ,Sp] = i€ a pjSj, and v S 



A total spin 



basis (TSB) \f m 



6 V satisfies the eigenvalue relations 

S%L) = J(J + l)fi>. S z \f m ) = m\%l), (1) 

and is associated with the direct sum decomposition of 
V into irreducible representations (irreps) of SU(2) [9(, 



(2) 



Here is a dj .-dimensional space that accounts for 

the degeneracy of the spin-j irrep and has basis |^) € 

where t = l,---,dj, whereas is a (2j + 1)- 
dimensional space that accommodates a spin-j irrep and 
has basis |^|) € where m is the projection of the 

spin in the z direction, m — —j, ■ ■ Each vector of 
the TSB factorizes into degeneracy and irrep parts as 
&> = OO, ^re Eq. © only determines ©. 

Bipartite decomposition.— A pure state |W) of a 
bipartite system with vector space A <g> B can always be 
expressed in terms of a TSB for A and a TSB for B as 
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When |W) is an SU(2) singlet, that is, invariant under 
transformations acting simultaneously on A and B, or 

(S [A] + S [B] ) 2 \V) =0, (S [A] +Si B] )\V) =0, (4) 
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then the symmetry materialises in constraints for the ten- 
sor of coefficients N, which splits into degeneracy and 
irrep parts according to [HI 
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where lo is completely determined in terms of Clebsch- 
Gordan coefficients (jxj2m\m2\ji32\ jm) [11], namely 
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FIG. 1: Diagramatic representation of tensors A and T of an 
MPS and tensors (77, w) and (X, C) of an SU(2) MPS. 



Eq. (JSJ) is quite sensible: it says that a coefficient N^^ml 
in Eq. ([3]) may be non-zero only if (i) j\ — ji (only the 
product of two spin j irreps can give rise to a spin 
irrep, that is, the singlet |\&)) and (ii) mi = — m,2, which 
guarantees that the z-component of the spin vanishes. In 
addition, Eq. (O embodies the essence of our strategy: 
to isolate the degrees of freedom that are not determined 
by the symmetry - in this case the degeneracy tensor 
T/ it;j . We now consider the singular value decomposition 



^ 2 = E(^)^M^') t 



of tensor Tl t for a fixed j, and define 
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By combining Eqs. ([5]), ([8]) and ([9]) we arrive to our 
canonical symmetric bipartite decomposition (CSBD) 



i*)=E(E^i< 1 )i< 1 )VE<iSl)ia 

which is related to the Schmidt decomposition 

I*) = ^a q |$L ai )|$L s1 ), 

ct 

by the identifications a — > (jtm), X a — * Vt^m an d 
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where some of the Schmidt coefficients A a are negative. 
More generally, a state |^^) of a bipartite system C®D 
can be expressed in terms of TSBs for C and D as [irj| 
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where tensor X relates degeneracy degrees of freedom 
and tensor C is given by the Clebsh-Gordan coefficients 



rum 

jl™. 1 j 2 m 2 



(jij2mim 2 \jij2;jm) 



(14) 



Matrix Product decomposition.— We now consider 
a chain of n quantum spins with spin s, represented by 
a ID lattice where each site, labelled by r (r = 1, . . . , n), 
carries a (2s + l)-dimensional irrep of SU(2). The coeffi- 
cients c mim2 ... mn of a state of the lattice, 



2s+l 



2s+l 



E-E 
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|!&>l&> ■••Ifel), (15) 



mi — 1 



where {|m)} is a basis for site r with Sz 
can be codified as an MPS P, Q , 



ai ■ ••a rl — \ 



r [l]mi A [l] r [2]m 2A [2] . , _y 



(16) 



Following the conventions of Q , here Xa ' are the Schmidt 
coefficients of according to the bipartition [1 • • -r] : 
[r + 1 • • • n] of the spin chain, while tensor rjTg™ relates 
Schmidt vectors for consecutive bipartitions, 

2s+l 



] ) = e r w i& ] >i*r"" w] >- 



m— 1 

When |^) is a singlet, that is 



E4 r] w = o, 



(17) 



(18) 



then Eqs. (JTOJ) and |(T3J) supersede Eqs. dTTJ> and (JT7J) 
and each tensor A and T in Eq. ([To]) decomposes into 
degeneracy and irrep parts, see Fig. (fT]), 
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Xjtj'v Cj, m , sm ,,, (20) 



where C 1 is related to Clebsch-Gordan coefficients C by 

< 5j'm'i"m" = (~1) 2J ( w m') ^j'™ " ( 2i ) 
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FIG. 2: The TEBD algorithm is based on updating the MPS 
when a gate U acts on two neighboring sites. This diagramm 
generalizes Fig. (3.i) in [l2T ] after the replacements A — » (rj, w) 
and T -» (X, C) of Eqs. ((IlJl-lgOl) for an SU(2) MPS. 



The SU(2) MPS is denned through Eqs. P])-([2Tj|l. In 
this representation, the constraints imposed by the sym- 
metry are used to our advantage. By splitting tensors A 
and r, we achieve two goals simultaneously. On the one 
hand, the resulting MPS is guaranteed, by construction, 
to be invariant under SU(2) transformations. That is, 
any algorithm based on this representation will preserve 
the symmetry exactly and permanently. On the other 
hand, all the degrees of freedom of \^f) are concentrated 
in smaller tensors r\ and X (tensors u> and C are speci- 
fied by the symmetry), and thus the SU(2) MPS is a more 
economical representation. If | • | denotes the number of 
coefficients of a tensor, then 

|A|=£(2j + l)d 3 - -> M = 5^>, (22) 

3 3 

|r| = (2 S +l)|A||A'| - \X\=J2Wi', (23) 

where A and A' are the tensors to the left and to the right 
of r, and where, following spin composition rules, the last 
sum is restricted to pairs (j, f) such that \j — j'\ < s. 

Simulation of time evolution.— Our next step is 
to generalize the TEBD algorithm 0] to the simulation 
of SU(2)-invariant time evolution. This reduces to ex- 
plaining how to update the SU(2) MPS when an SU(2)- 
invariant gate U acts between contiguous sites, see Fig. 
([2]). The update is achieved by following steps analo- 
gous to those of the regular TEBD algorithm, see Fig. 
(3) of |l2|, involving tensor multiplications and one sin- 
gular value decomposition (SVD), Fig. ([3]). However, 
all these manipulations involve now smaller tensors, and 
only tensors X and r\ of the SU(2) MPS need to be up- 
dated. This results in a substantial reduction of compu- 
tational space and time, and thus an increase in per- 
formance. For instace, the SVD of O in Fig. (3) of 
where |6| s» (2s + 1) 2 |A| 2 , is now replaced with the 
SVD of Sljtjt' ( see Fig. ([3])) for each value of j, where 
l%«t'l = (J2^>j-2s d 3') 2 - Tne cost c svd(A) of comput- 
ing the SVD of a matrix A grows roughly as |A| 3 / 2 and 
is the most expensive manipulation of the TEBD algo- 




FIG. 3: Key step of the TEBD algorithm for an SU(2) MPS, 
analogous to Figs. (3.i)-(3.m) in \12\ for a regular MPS. — 
Once U has been applied on two spins, additional tensors 
Vxi, Vci implement a unitary transformation required to re- 
absorb these spins into blocks and obtain an updated repre- 
sentation for the bipartition [1 • • • r] : [r + 1 • • • n]. Then, for 
each fixed value of the j indices (discontinuous lines), the 77's, 
X's and Vx's are multiplied together and the result, with a 
weight coming from the product of the ui's, C's, Va'a and 
U [that can be pre-computed because none of these tensors 
depend on |^/}], is added together to give rise to tensor Q. 
A singular value decomposition of Cljtjt' f° r each value of j 
ensues, see Eq. (|25[1 . and minor rearrengements finally lead 
to updated tensors X' [r \ r/' lr] and X' [r+1] . 



rithm. We obtain the following comparative costs: 

c svd (&) ~ [(2* + l)5}(2j + l)d 3 -]j , (24) 

(25) 



j+2s 



Example.— For illustrative purposes, we consider a 
quantum spin chain with s = 1/2 and with Hamiltonian, 

H = Y / (S l I ] ^ r+1] + 4 rI 4 r+1] + SP St +1] ) , (26) 

r 

that is, the spin-1/2 antiferromagnetic Heisenberg model, 
which is SU(2) invariant and quantum critical at zero 
temperature. We have computed an SU(2) MPS approx- 
imation to the ground state of H , in the limit n — > 00 of 
an infinite chain, by simulating imaginary-time evolution 
[l2| starting from a state made of nearest-neighbor sin- 
glets (\ [ l] 2 )tt}\) -|-i/ 2 )lK 1] ))/^- With the constraint 
J2j dj — 600, we have obtained that the following irreps 
j, with degeneracies dj, contribute to the odd and even 
bipartitions of the resulting state, 
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Eqs. (|22]) - ([25f show substantial computational gains, 
|r| 10 7 „ c svd (Q) 9xl0 10 



\x\ 



2 x 10 



5 



50, 



3 x 10 8 



300, (27) 



(O) 

that is, with a regular MPS, storing the same state would 
require about 50 times more computer memory, while 
performing each SVD would be about 300 times slower. 
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FIG. 4: (Up) Errors in the two-point correlator C2(r) for 
1 < r < 7, when using an SU(2) MPS of different sizes 
Xi G {350, 700, 1110, 1450, 1800, 2200}. Here x is roughly |A| 
in Eq. (|22p . that is, the rank of an equivalent (regular) MPS. 
The lowest line, X6 ~ 2200, shows the errors in the data pre- 
sented in the table. (Down) Numerical results for 62(7") for 
up to r — 20,000 sites, for different sizes Xj> together with 
corresponding errors e^. 



We have computed the two-point correlators (r) = 
(S l " ] S l " ] ) and C 2 v (r) = (S [ ' ] S [ I +1] ) 0, and the average 
C 2 (r) = (C 2 A (r) + C 2 v (r))/2 Q. For small r they read: 
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C£(r) 


C 2 V (r) 


Ca(r) 


1 


-0.14800224748 


-0.14742920605 


-0.147715726[7] 


2 


0.06067976982 


0.06067976991 


0.060679769[9] 
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-0.05037860908 


-0.05011864581 


-0.050248627[4] 
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0.03465277614 


0.03465277645 


0.034652776[3] 
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-0.0309785296 


-0.0308021901 


-0.03089036[0] 
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0.024446726 


0.024446726 


0.0244467[26] 


7 


-0.022565932 


-0.022430482 


-0.0224982[1] 



where, for C?2(r), the square brakets show the first digits 
that differ from the exact solution [Tsl ] , from which we re- 
cover e.g. 9 significant digits for r = 1. An expression for 
the correlator C 2 (2) is also known for large r 17J. There, 
for r w 4, 000, 10, 000 and 13, 000, our results approxi- 
mate the asymptotical solution with an error of 1%, 5% 
and 10% respectively, see Fig. (|4|). For comparison, with 
a regular MPS and similar computational resources, we 
lose three digits of precision for r = 1, whereas a 10% er- 
ror is already achieved for r w 500 instead of r w 13, 000. 

Final Remarks.— The above test with a critical spin- 
1/2 chain unambiguously demonstrates the superiority 
of the SU(2) MPS and TEBD with respect to their non- 
symmetric versions. Most promissingly, these techniques 
can now be used to address systems that remain other- 
wise largely unaccessible to numerical analysis due to a 
large dimension of the local Hilbert space. These include 



a chain made of large spins, say s — 4, or a spin ladder 
with several legs. We regard the latter as a chain with 
several spins per site, where each site decomposes into 
SU(2) irreps as in Eq. © [3]. 

In addition, the SU(2) MPS is not restricted to the 
representation of SU(2) singlets. On the one hand, it can 
be used to represent any SU(2) invariant mixed state p 
of the chain, which decomposes as (see Eq. (|5J)) 



= 0Pj<8ii 



2j+l 



(28) 



This is achieved by attaching, to the end of the chain, 
an environment E that duplicates the subspace V of the 
chain on which p is supported, and by considering a sin- 
glet purification \$>J E ), where p = tr E \^ E )(^ E \. We 
first build an SU(2) MPS for the purification and then we 
trace out E. The resulting structure is a matrix product 
representation that retains the advantages of the SU(2) 
MPS. In particular, notice that when p corresponds to 
one single irrep j, 
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P = 



2j 



- y \ [v] )i lv] \ 

/ j \jm/\jm\ 
m=-j 



(29) 



then the environment is a site with a spin j, and the 
chain together with the environment is just an extended 
spin chain, with the purification being of the form 



1 



y ' \jm) 
m=-j 



(30) 



On the other hand, the SU(2) MPS can also be modified 
to represent any pure state |^|) of the chain with well 
defined j and m. To see this, we first consider a mixed 



state p as in Eq. (f2"9"| . that is, a symmetrization of 
and then a purification \^ p } for p as in Eq. (|29[) . for 
which we can build an SU(2) MPS. Finally, we recall that 
|W) = which leads to a simple, SU(2) MPS- 

like representation for in terms of the SU(2) MPS 
for the purification \^ p ). The time-evolution simulation 
techniques described in this paper can be applied to the 
above generalized representations. 

Near the completion of this paper, we became aware 
of related results by I. McCulloch derived independently 
in the context of DMRG [ll|. 

The authors thank S. Lukyanov for helpful communi- 
cations. G.V. acknowledges support from the Australian 
Research Council through a Federation Fellowship. 
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For semi-integer spins, e.g. s = i , the CSBD of Eq. (Tl0|) 
for a partition [1 • • ■ r] : [r+ 1 • • • n] of the chain has only 
integer (semi-integer) values of j when r is even (odd). 
CJ^l) (w =A, V) are computed by contracting a small 
tensor network involving the tensors (r/, lj) and (X, C) for 
two contiguous sites. For r > 1, C_"(r) is computed in the 
same way, but first we simulate r — 1 (SU(2) invariant) 
swap gates that bring the two relevant sites together. 
We find, numerically, that the effect of [l_| persists in the 
thermodynamic limit (open boundary conditions). As a 
result, the ground state is invariant under shifts by two 
(but not one) lattice sites, and C% (r) 7^ C% (r) for odd r. 
See J. Sato, M. Shiroishi and M. Takahashi, Nucl. Phys. 
B 729, 441 (2005), and references therein. 
We use c = 1.05 in Eq. (5.25) of S. Lukyanov, V. Terras, 
Nucl. Phys. B654 (2003) 323-356, |hep-th/0206093"l 
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